1. Field of the Invention
This invention relates generally to the field of geophysical prospecting. More particularly, the invention relates to the field of seismic data processing. Specifically, the invention is a method for downward extrapolation of pre-stack seismic data for seismic migration.
2. Description of the Related Art
Seismic migration is the process of constructing the reflector surfaces defining the subterranean earth layers of interest from the recorded seismic data. Thus, the process of migration provides an image of the earth in depth or time. It is intended to account for both positioning (kinematic) and amplitude (dynamic) effects associated with the transmission and reflection of seismic energy from seismic sources to seismic receivers. Although vertical variations in velocity are the most common, lateral variations are also encountered. These velocity variations arise from several causes, including differential compaction of the earth layers, uplift on the flanks of salt diapers, and variation in depositional dynamics between proximal (shaly) and distal (sandy to carbonaceous) shelf locations.
It has been recognized for many years in the geophysical processing industry that seismic migration should be performed pre-stack and in the depth domain, rather than post-stack, in order to obtain optimal, stacked images in areas with structural complexity. In addition, pre-stack depth migration offers the advantage of optimally preparing the data for subsequent AVO (Amplitude Versus Offset) analysis. Pre-stack depth migration has traditionally been performed through the application of Kirchhoff methods. However, because of recent advances in computing hardware and improvements in the design of efficient extrapolators, methods that are based on solutions of the one-way wave equation are becoming increasingly popular.
It has been well established in the literature that wave equation-based methods offer the kinematics' advantage of implicitly including multipathing, and can thus more accurately delineate structures underlying a complex overburden. However, there has been considerably less discussion of the dynamic advantages that might (or should) be realized through wave equation-based methods. This is not surprising, since Kirchhoff pre-stack depth migration has undergone a much longer evolution than its wave equation-based counterparts. Part of this evolution has been the development of various factors or strategies to account for geometrical divergence and spatial irregularities in acquisition.
One method of prestack depth migration individually treats common shot gathers, or shot records. This method is thus called “shot record” or “shot profile” migration. Shot record migration is capable of providing excellent imagining accuracy and excellent amplitude preservation. These advantages have made shot record migration one of the more popular methods of wave theory migration.
An early implementation of shot record migration is described in Reshef, Moshe and Kosloff, Dan, “Migration of common shot gathers”, Geophysics, Vol. 51, No. 2 (February, 1986), pp. 324–331.
Seismic migration comprises two steps: (1) wave extrapolation and (2) imaging. The present application deals with the first step of downward extrapolation in seismic migration. A corresponding method of imaging in seismic migration is described in co-pending U.S. patent application, “Method for Imaging of Prestack Seismic Data”, by the same inventor and assigned to the same assignee as the present patent application.
Downward wave extrapolation results in a wavefield that is an approximation to the one that would have been recorded if both sources and receivers had been located at a depth z. One way in which seismic migration methods differ is in their method of implementation of the wave extrapolation step, due to differences in solving the wave equation.
Gazdag extrapolation is described in Gazdag, Jenö, “Wave equation migration with the phase-shift method”, Geophysics, Vol. 43, No. 7 (December, 1978), pp. 1342–1351. Gazdag (1978) discloses a migration method for zero-offset seismic data, based on numerical solution of the wave equation using Fourier transforms. In the Gazdag extrapolation method, the source and receiver positions are extrapolated downwards by means of a phase shift, a rotation of the phase angle of the Fourier coefficients of the seismic data in the frequency domain. The computations are equivalent to multiplication by a complex number of unit modulus, assuring unconditional stability. Additionally, the method handles high propagation (dip) angles. However, the Gazdag phase-shift method only works well for laterally uniform velocity.
The Gazdag phase-shift method is extended to media with lateral velocity variation in Gazdag, Jenö, and Sguazzero, Piero, “Migration of seismic data by phase shift plus interpolation”, Geophysics, Vol. 49, No. 2 (February, 1984), pp. 124–131. Gazdag and Sguazzero (1984) discloses a two-step extrapolation method typically reffered to as phase shift plus interpolation. In the first step, the wavefield is extrapolated by the Gazdag phase-shift method using a plurality of laterally uniform velocity fields, generating reference wavefields. In the second step, the actual wavefield is computed by interpolation from the reference wavefields. As for the Gazdag phase-shift method, the Gazdag phase-shift plus interpolation method is unconditionally stable and can accurately handle very steep dips. The overall accuracy of migration increases with the number of reference velocities used, but so does the computational cost. Areas with very large lateral velocity variations can require up to five reference wavefields, each entailing a pair of spatial Fourier transforms. Thus, the Gazdag phase shift plus interpolation method can be very expensive.
Split-step Fourier extrapolation (also called the phase-screen method) is a modification of the Gazdag phase shift method to accommodate lateral velocity variations in each migration interval. Split-step Fourier extrapolation is described in Stoffa, P. L., Fokkema, J. T., de Luna Freire, R. M., and Kessinger, W. P., “Split-step Fourier migration”, Geophysics, Vol. 55, No. 4 (April, 1990), pp. 410–421. Stoffa et al. (1990) discloses a two phase-shift method to handle lateral velocity variation by first defining a reference slowness (reciprocal of velocity) as the mean slowness in the migration interval and secondly defining a source contribution by a spatially-varying perturbation or correction term. The first phase shift is identical to the constant velocity Gazdag phase shift method. The mean slowness defines a reference vertical wavenumber that is used to downward continue the data in the frequency-wavenumber domain. The second phase shift is a correction term providing a time shift based on the differences between the actual and reference slownesses at each spatial position. The two phase shifts are repeated for each migration interval. Split-step Fourier extrapolation works well for a broad range of incidence angles, provided that the relative velocity variation is small. It can also accommodate large velocity variations, but only for a limited range of angles.
Error analysis for the split-step Fourier method is described in Huang, Lian-Jie and Fehler, Michael, “Accuracy analysis of the split-step Fourier propagator: Implications for Seismic modeling and Migration”, Bulletin of the Seismological Society of America, Vol. 88, No. 1 (February, 1998), pp. 18–29. The split-step Fourier method expands the square root operator and splits the exponential operator in the one-way scalar wave equation. Huang and Fehler (1998) show that these approximations give an error under 5% for propagation angles less than about 45° and for relative velocity variations of less than about 10%. For small relative velocity variations, the split-step Fourier method is accurate for propagation angles up to about 60°.
In laterally homogeneous media, exact wave extrapolation can be achieved by applying single phase shift, which is a solution of the one-way wave equation in the frequency-wavenumber domain. The wave equation in this domain is sometimes referred to as the “square root equation.” In laterally heterogeneous media, however, wave extrapolation is achieved by applying a rational approximation of the square root equation. The approximation is typically carried out by finite difference methods. The two methods of rationalizing the square root equation are by Taylor series and by continued fractions. The Taylor series method is called an explicit scheme while the continued fraction method is called an implicit scheme. Implicit finite difference extrapolation methods are particularly well known in the art.
Implicit finite difference extrapolation methods are described, for example, in the two articles: Lee, Myung W. and Suh, Sang Y., “Optimization of one-way wave equations”, Geophysics, Vol. 50, No. 10 (October, 1985), pp. 1634–1637 and Li, Zhiming, “Compensating finite-difference errors in 3-D migration and modeling”, Geophysics, Vol. 56, No. 10 (October, 1991), pp. 1650–1660. Lee and Suh (1985) apply least squares to optimize coefficients of second order approximations to the square root equation. Li (1991) compensates for crossline and inline splitting errors in 3D migration by a phase-shift compensation filter applied ever few steps of finite difference extrapolation. The phase-shift compensation filter is based on the Gazdag phase shift or phase shift plus interpolation methods.
A hybrid migration method is described in Ristow, Dietrich and Ruhl, Thomas, “Fourier finite-difference migration”, Geophysics, Vol. 59, No. 12 (December, 1994), pp. 1882–1893. Ristow and Ruhl (1994) combines a Gazdag phase-shift extrapolation in the frequency-wavenumber domain and a subsequent finite difference extrapolation in the frequency-space domain. This method combines the split-step Fourier method with the finite difference method. The finite difference step corrects for the effect of large lateral velocity perturbations on the diffraction term. In this respect, the finite difference complements the preliminary perturbation in the spatial domain for thin lens propagation.
If the lateral velocity variations are sufficiently small, then the preliminary, constant-velocity phase shift accounts for the majority of the moveout in the diffraction term. Thus, the finite difference computation may not require any phase correction. Under these circumstances, for each wavefield, there is only one Gazdag phase shift in the (ω, k) domain, followed by one thin lens phase shift and one finite difference computation in the (ω, x) domain. However, if the lateral velocity variations are sufficiently large, then one or more phase corrections may be required after finite differencing, each of which entails forward and inverse spatial Fourier transforms and a phase shift. The cost of these phase shifts is in addition to the cost of the preliminary Gazdag phase shift.
A recent implementation of this hybrid migration method is discussed in Sun, James, Notfors, Carl, Gray, Sam and Zhang, Yu, “3D Pre-stack common shot depth migration: A structural adaptive implementation”, SEG Int'l Exp. & Ann. Mtg., San Antonio, Tex., Sep. 9–14, 2001, Expanded Abstracts, pp. 1017–1020. Sun et al. (2001) describe a procedure for deciding between the Gazdag and finite difference methods for each migration interval. They also develop their own sets of coefficients for the finite differencing. Still, the computation method for the finite difference coefficients appears to be somewhat arbitrary and difficult. In addition, no authors have provided an algorithm that quantitatively identifies the structural conditions under which the trailing finite difference computation is required.
Thus, there exist extrapolation algorithms that are appropriate for homogeneous, weakly heterogeneous, and strongly heterogeneous media, as manifested by degree of lateral velocity variation. Not surprisingly, the more accurate methods are also the more costly. It would thus be desirable to be able to identify the most cost-effective method for a particular extrapolation step. Thus, a need exists for a hybrid migration method that chooses the most efficient extrapolator for each depth step by an automated error analysis based on the complexity of the underlying velocity model.